The Journal of Chemical Physics

Synapsin-1 is a multidomain neuronal protein containing extensive intrinsically disordered regions and is a key component of synaptic-vesicle condensates. Direct residue-level simulation of the collective organization of thousands of synapsin molecules remains computationally demanding. Here, we develop a coarse-grained description that connects residue-level CALVADOS 3 simulations to a one-particle-per-protein model. A potential of mean force between two synapsin molecules is obtained by umbrella sampling and represented by an isotropic effective interaction containing a short-range attractive region and a weak outer repulsive contribution. We compare two treatments of this interaction that differ only in the retention of the outer tail. Langevin dynamics simulations of effective proteins show aggregation upon cooling and compression in both models, but with markedly different collective organization. The shorter-ranged model progressively coarsens toward a single dense domain, whereas retaining the outer repulsive contribution favors the persistence of multiple mesoscale aggregates. The two models also display distinct relationships between aggregate size and particle mobility at low temperature. These results show that weak features of an effective protein-protein interaction can have pronounced consequences for collective synapsin organization at mesoscopic scales. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=95 SRC="FIGDIR/small/745062v1_ufig1.gif" ALT="Figure 1"> View larger version (23K): org.highwire.dtl.DTLVardef@854288org.highwire.dtl.DTLVardef@d2fc52org.highwire.dtl.DTLVardef@1b37d5aorg.highwire.dtl.DTLVardef@eac583_HPS_FORMAT_FIGEXP M_FIG C_FIG

The process of drug unbinding is of immense importance in the field of biophysics and therapeutics. The behavior of these systems is greatly influenced by their thermodynamic and kinetic properties. Therefore, it is crucial to accurately estimate the ligand binding free energies and rate of ligand dissociation, yet these processes are often governed by rare event transitions that lie beyond the reach of standard brute-force molecular dynamics simulations. While enhanced sampling simulations offer a solution, their efficacy is strictly contingent upon the selection of appropriate collective variables (CVs) which is non-trivial for complex systems like protein-ligand complexes. In this study, we present a method to derive optimized CV from transition state region (TS) via an interpretable machine learning (ML) model, Elastic Net. By employing some physically intuitive order parameters, the derived optimized CV from the TS-region greatly accelerate ligand binding-unbinding transitions and achieves rapid free energy surface (FES) convergence across diverse systems including buried and solvent exposed active sites such as Trpsin-benzamidine complex, host-guest systems and sodium epoxidase etc. Intriguingly significant contribution of the ligand hydration is found in the optimized CV which depicts crucial role of solvent in driving ligand binding-unbinding transitions. The estimated binding free energies for different protein-ligand complexes match quite well with experiments, while maintaining a low computational cost. The derived optimized CV is also used to calculate the ligand residence times across different systems and calculated residence times are within the experimental range for all systems, again with very little computational costs. Moreover, we show that the optimized CV constructed from TS region via an interpretable ML model is transferable across diverse systems, offering a robust and scalable framework for drug discovery and investigation of complex biomolecular recognition.

Gaussian accelerated molecular dynamics (GaMD) enhances conformational sampling by adding a smooth boost potential without requiring predefined collective variables, but an engine-integrated implementation has not been available in GROMACS. Here, we implement total-, dihedral-, and dual-boost GaMD in GROMACS 2025.4, including staged energy-statistics collection, GPU-based bias evaluation and force scaling, restart support, and outputs required for cumulant-based free-energy reweighting. The implementation was evaluated using four benchmark systems spanning conformational free energies, protein folding, and ligand recognition. For alanine dipeptide, a reweighted 100 ns GaMD trajectory recovered the major free-energy basins and rotational barriers in overall agreement with a 1000 ns conventional MD simulation. For chignolin and TC5b, all three independent trajectories for each system sampled native-like folded states from extended conformations within 300 ns and 1 s, respectively; the best TC5b structure had a minimum backbone RMSD of 0.03 nm from the experimental structure. In the benzene-T4 lysozyme system, two of five independent 500 ns trajectories captured both ligand binding and dissociation, yielding a bound pose with a minimum ligand RMSD of 0.06 nm from the crystal structure. Across all four systems, the boost-potential distributions were approximately Gaussian, and second-order cumulant reweighting resolved the expected conformational and binding free-energy basins. These results demonstrate that GROMACS-GaMD provides a practical, GPU-enabled, collective-variable-free enhanced-sampling framework for biomolecular free-energy calculations, protein folding, and ligand-binding studies.

Recent experiments show that protein condensates sitting on opposite surfaces of a flat lipid membrane move together and prefer to overlap, even though they cannot touch each other. This points to an indirect, membrane-mediated interaction. Two mechanisms could be responsible: a curvature-induced interaction, which is energetic in origin, and a fluctuation-induced interaction, which is entropic. Here we study both with coarse-grained molecular dynamics simulations, using Cookes implicit-solvent lipid model together with a generic bead-spring polymer model for the condensate. We compute the potential of mean force between two condensates across the membrane. For condensates of the same size, full overlap is unfavorable, and the pair instead settles into a partially overlapping state that bends the membrane into an S-like shape. When the two condensates differ strongly in size, full overlap becomes favorable. We explain this with a simple geometric picture. The condensate wets the membrane as a thin film and imposes curvature only along its rim, while membrane tension flattens the membrane under its interior. The resulting ring of curvature can trap a smaller condensate on the opposite side. We also compare the bending undulations and the effective bending modulus of a bare membrane, a membrane with one condensate, and a membrane with condensates on both sides. A wetting condensate suppresses the undulation modes and stiffens the membrane, but whether this makes overlap entropically favorable remains inconclusive. Our results indicate that the coupling is driven mainly by curvature, and that it depends on the wetting mechanism and on the membrane tension.

Polyglutamine (polyQ) aggregation plays a central role in several neurodegenerative diseases, including Huntington's disease. DNAJB6b, a molecular chaperone involved in protein quality control, is known to efficiently suppress polyQ aggregation, but its anti-aggregation mechanism remains unclear. In this work we investigate the interaction between DNAJB6b and the polyQ region (Q48) of mutant Huntingtin Exon 1 (mHttEx1) using a custom-built coarse-grained molecular dynamics model. The model incorporates a two-bead-per-amino-acid representation with hydrogen bonding (termed 2BPA-HB), and is calibrated against all-atom molecular dynamics data in terms of geometry, hydrophobicity, and hydrogen bonding. The model reproduces the tertiary structure of DNAJB6b and its interactions with Q48, and reveals an inverse correlation between DNAJB6b concentration and Q48 aggregation propensity. Our simulations show that DNAJB6b co-condensates with polyQ molecules, thereby shielding the polyQ from forming the intermolecular hydrogen bonds necessary for amyloid formation. The 2BPA-HB CGMD model en- ables efficient exploration of DNAJB6b conformations, supporting future studies of chaperone-mediated aggregation suppression and therapeutic development.

In the active loop extrusion model, the cohesin protein complex creates chromatin loops in eukaryotic cells. Extrusion maintains topologically associated domains (TADs), which are contiguous segments of chromatin that preferentially colocalize in space and are typically bounded by CTCF proteins that pause cohesin translocation. Here, we model active loop extrusion with hybrid molecular dynamics - Monte Carlo simulations in entangled flexible linear polymer melts. Intra-chain contact probabilities of polymers with active loop extrusion are enhanced compared to their equilibrium, passive counterparts. Extrusion causes the size of chain segments to be much smaller than in passive melts. While the overlap parameter in passive melts without extrusion monotonically increases with segment length, it is nonmonotonic in active melts and on the order of unity within the parameters of this study. Active loop extrusion suppresses contacts between TADs in favor of intra-TAD contacts. Reduction of overlaps between chain segments dilutes entanglements in active melts. Depending on parameters, active extrusion without TADs may induce more compact conformations than with TADs, due in part to fractal loopy globule-like dynamics. This work suggests that active loop extrusion reduces overlaps between TADs, contributing to effective gene regulation by cis-regulatory elements.

Interleukin-13 (IL-13) is an immunomodulatory cell signaling cytokine that has been implicated in neurodegenerative disease and chronic inflammation. IL-13 binds to its low and high affinity receptors, IL-13 receptor 1 (IL-13R1) and IL-13 receptor 2 (IL-13R2), respectively, with residence times that vary accordingly. As the binding kinetics of the cytokine-receptor complexes influence cellular responses, we employed the molecular dynamics (MD) simulation-based{tau} -random acceleration molecular dynamics method ({tau}RAMD) to compute relative residence times for wild-type (WT) IL-13 and 19 IL-13 mutants to the two receptors. Comparison with experimental kinetic data shows that the{tau} RAMD computations capture the trends in residence times. Analysis of simulated dissociation trajectories of the cytokine-receptor complexes reveals two distinct dissociation pathways of IL-13 from each of the receptors. This study thus pinpoints key determinants of the interaction of IL-13 with its receptors which could be targeted for therapeutic design. Statement of SignificanceCytokines are regulatory proteins that bind to cell surface receptors and thereby send signals to the cellular interior. Interleukin-13 (IL-13) is a cytokine that has a low and a high affinity receptor. It has important physiological roles, and its deregulation is involved in diseases such as atopic dermatitis and asthma. We employed a molecular dynamics simulation-based method to compute the effects of changes in the sequence of IL-13 on the lifetimes of complexes of IL-13 and its receptors. Comparison with experiments supports the validity of the computational approach and analysis of the simulations reveals two distinct ways in which IL-13 dissociates from each receptor. These results thus provide a map for targeting IL-13 - receptor interactions for the design of therapeutics.

Neutron scattering is a powerful technique to investigate atomic structures and molecular dynamics of proteins at the nano-scale. When it comes to dynamics, incoherent and coherent scattering respectively provide information on the single and collective dynamics of nuclei. In proteins, hydrogen has the highest incoherent cross-section, and it is common practice to overlook the contribution of coherent terms stemming from all nuclei. However, the fast collective dynamics of heavier nuclei could also be studied if coherent scattering and incoherent scattering were experimentally separated. The recent advent of polarized neutron spectroscopy with sufficient flux and energy resolution has made it possible, and opens new perspectives to investigate the relative importance of coherent scattering and the information it provides on biological samples. The present study reports on the use of polarized quasi-elastic neutron scattering (QENS) and the application of a minimalistic model adapted to both individual and collective dynamics. Using a perdeuterated green fluorescent protein as a model globular protein, the study provides an interpretation of the dynamical parameters obtained with QENS, and a comparative study of the Elastic Coherent and Incoherent Scattering Factor. Based on both experiments and calculations, we discuss the relative importance of distinct and self components of coherent scattering, which is often wrongly assumed to be representative of collective dynamics only. The results highlight the current impediments rendering complicated a straightforward analysis of fast collective dynamics in hydrated protein samples.

Calcium oxalates are naturally occurring minerals, with the monohydrate form CaC2O4xH2O (COM) being the most stable. COM contains two crystallographically different water sites. We investigate the details of water internal dynamics in the high and low temperature phases of COM using 17O central transition solid-state NMR line shapes, as well as laboratory and rotating frame relaxation rates. The measurements were performed either under static or magic angle spinning conditions and in a wide temperature range from 343 to 180 K. The combination of all measurements allows for precise constraints on motional mechanisms, rate constants, and amplitudes of motions. The high temperature phase is dominated by large-angle fluctuations with an amplitude of about 100 degrees, identical in both sites. During the phase transition between 323 to 300 K, these large-angle jumps freeze out in one of the water sites, while remaining active in the other. In the low temperature phase from 280 to 180 K, small-angle fluctuations of 2-8 degrees in amplitude dominate the relaxation. Transverse relaxation rates also point to the existence of a very slow collective rocking motion down to about 220-200 K.

Contact inhibition of proliferation (CIP) enables epithelial tissues to self-regulate growth and maintain tissue homeostasis. However, how cell-level mechanical contact, tissue-scale structural order, and proliferation kinetics interplay remains a fundamental open question in living matter physics. Here, we present a particle-based model of a confluent epithelial monolayer governed by overdamped dynamics, where individual cells interact via a two-dimensional hard core- soft shoulder potential. By comparing structural evolution during quasistatic densification with previously reported experimental division kinetics, we find that the dynamics of proliferation arrest mimics the onset of direct steric contacts between the hard cores of the shell. Identifying hard core contacts as the physical driver of CIP, we couple our mechanical model with a stochastic Monte Carlo division scheme in which the instantaneous division rate decreases to zero from an intrinsic value as the number of hard core contact increases to six from zero. We demonstrate that for high intrinsic division rates, the cellular densification outpaces mechanical relaxation. This kinetic mismatch drives premature hard-core contact formation, shifts the onset of jamming and contact inhibition to lower packing fractions, and induces increasingly disordered transient configurations before the tissue universally converges to a hexagonal close-packed limit. Our model's predicted division kinetics and structural order evolution are consistent with epithelial monolayer experiments, both reported and our own. This minimal physical framework links single-cell steric contact mechanics directly to tissue-scale growth regulation and structural evolution.

The dynamics of gene regulatory networks are governed by intrinsic noise, stemming from the random nature of biochemical reactions, and by extrinsic noise, arising from fluctuations in cellular components and environmental conditions. Together, these sources can compromise the reliability of predictive computational models if not properly accounted for, and capturing both effects within a single framework remains a non-trivial task in computational biology. In this work, we propose an uncertainty quantification framework that addresses these two contributions jointly: intrinsic stochasticity is described through a partial integro-differential equation (PIDE) for the protein probability density function, whereas extrinsic noise is represented as parametric uncertainty in the kinetic parameters. The propagation of the uncertainty is carried out via an intrusive polynomial chaos expansion (PCE), in which the PCE coefficients are obtained from a stochastic Galerkin projection of the PIDE, yielding a coupled deterministic system that is solved with standard numerical methods. We illustrate the approach on a positive autoregulatory gene network with one and two uncertain kinetic parameters. The proposed approach accurately reproduces the mean, variance, and full protein probability density function, including the bimodal distributions, at a substantially lower computational cost.

A fundamental biological process of transcription occurs in the cell nucleus, which is a complex medium that also contains multiple heterogeneous structures known as biomolecular condensates. Interestingly, some of these condensates contain mRNA molecules in addition to proteins, suggesting an important cellular role in transcription that is not yet well understood. In this work, we develop a minimal theoretical framework for quantitative investigation of the role of reversible mRNA condensation in transcription. Our discrete-state stochastic approach accounts for the most relevant processes, allowing us to explicitly evaluate the properties of the system and clarify the effects of condensation. Analytical calculations supported by computer simulations suggest that reversible mRNA condensation influences the transcription processes by maintaining a constant level of free mRNA in the nucleoplasm while lowering the degree of stochastic noise and increasing the robustness against external perturbations. Physicochemical arguments are presented to explain these observations. The proposed theoretical framework elucidates important microscopic aspects of transcription, providing a convenient quantitative tool for investigating complex biological phenomena.

Biomolecular condensates function as membraneless compartments, and some protein condensates can selectively concentrate single-stranded nucleic acids while excluding double-stranded nucleic acids. Understanding how nucleic acid structure affects partitioning into condensates has important implications for nucleic acid activity and function within condensates. Here, we present a set of coarse-grained two-bead-per-nucleotide models for simulations of double-stranded RNA and DNA in the CALVADOS framework. Our models separately represent the backbone and base, and maintain the helical structures using an elastic network potential tuned to capture chain stiffness. For dsRNA, the base stickiness was tuned using experimental data on differential partitioning of single- and double-stranded RNA into Ddx4N1 condensates in order to account for reduced base accessibility upon duplex formation. This RNA structural selectivity varied with the balance of electrostatic and non-electrostatic interactions, as revealed by simulations of condensates of the CAPRIN1 disordered region at varying ionic concentrations and with an R-to-K sequence variant. Finally, we developed parameters for double-stranded DNA using a similar approach. We envision that the CALVADOS models for double-stranded RNA and DNA will be useful for studying co-condensates of proteins and structured nucleic acids.

High-dimensional population balance equations (PBEs) provide a natural framework for modeling heterogeneous cell populations, but their direct numerical solution becomes computationally prohibitive when the internal state space contains many molecular variables. We propose a hybrid mechanistic-machine learning framework for reducing and simulating PBEs defined over high-dimensional intracellular coordinates. The cell population is described by a number density n(x, t), where x [isin] [R]N represents gene and protein states associated with macrophage activation. A dynamics-preserving autoencoder maps this state space to a low-dimensional latent coordinate z [isin] [R]d, with d << N, while retaining key qualitative features of the underlying gene regulatory network, including attractor structure and multistability. Mechanistic information from the original regulatory dynamics is used to construct interpretable drift and diffusion terms for the reduced latent-space PBE. The reduced PBE is solved using a stochastic Lagrangian particle representation, in which particles evolve according to stochastic differential equations (SDEs) corresponding to the latent drift and diffusion fields. The resulting latent-space solution is subsequently decoded and propagated back into the original state space to recover physically interpretable cellular dynamics. We demonstrate the framework on macrophage polarization under cytokine-dependent regulation, including gene knockout perturbations. Overall, the proposed framework provides a computationally tractable and mechanistically interpretable route for integrating single-cell genomic data with population balance models of cell-state dynamics.

This present work shows that E-hook fragments possess functional structure differences governed by electrostatic interactions and sequence composition. The acidic C-terminal tails of tubulin, known as E-hooks, play a central role in regulating interactions between microtubules and motor proteins, microtubule-associated proteins, and enzymatic modifiers. Despite their functional importance, the intrinsic structural properties of these peptide segments remain poorly characterized due to their intrinsically disordered nature. In this work, we present quantum-mechanically optimized structures of hexamer peptides derived from {beta}-tubulin E-hook sequences. Density functional theory calculations were used to optimize peptide geometries using progressively larger basis sets. From the optimized geometries we calculated theoretical Raman spectra, Ramachandran backbone dihedral distributions, and measured radii of gyration to resolve composition dependent structural tendencies. The combined Raman and conformational analyses provide a systematic computational approach for comparing simulated and experimental Raman spectra of tubulin E-hooks and other intrinsically disordered proteins and offer insight into how E-hooks contribute to the recognition mechanisms underlying the tubulin code.

Biomolecular condensates formed by intrinsically disordered proteins require molecular models that accurately describe proteins in both dilute solution and condensed phases. Explicit-solvent coarse-grained models offer an attractive balance between chemical resolution and computational efficiency. Yet, it remains unclear whether improving dilute-state properties is sufficient to obtain an accurate description of condensates. Here, we address this question by introducing minimal modifications to the Martini 3 force field that combine recent advances in bonded interactions with refined protein-water interactions and strengthened glycine self-interactions, while preserving the underlying chemical transferability of the model. The resulting model substantially improves the description of single-chain conformations across a diverse benchmark of disordered proteins. We then investigate phase separation of the well-characterized low-complexity domain of heterogeneous nuclear ribonucleoprotein A1 and its sequence variants. The model reproduces several key physicochemical properties of biomolecular condensates, including chain expansion in the dense phase, sequence-dependent intermolecular contacts, protein diffusion and its relation to single-chain dimensions, and hydration, while revealing quantitative limitations in condensate density, phase equilibria, and ion partitioning. Our results show that improving dilute-state behaviour translates into a better description of condensed-phase properties, including condensate density, but is not sufficient to quantitatively reproduce the equilibrium between the dilute and dense phases.

N-linked glycosylation is known to modulate the catalytic function of human carbonic anhydrase (HCA) IX, yet its influence on the underlying free-energy landscape remains largely unexplored. In the present work, we combine extensive all-atom molecular dynamics simulations with kinetic transition network analysis to investigate the effect of glycosylation on the conformational organization of the catalytic domain of HCA IX in both monomeric and dimeric forms. The multidimensional conformational space is discretized into distinct free energy minima using the distribution of reciprocal interatomic distances (DRID), and the effective barriers separating them are estimated using the max flow-min cut formalism. The corresponding free energy landscapes are visualized in terms of disconnectivity graphs, which provide a faithful representation of underlying kinetics. Minimum free energy paths, mean first passage times, as well as frustration metrics are computed to further quantify the effect of glycosylation on landscape topography. Unglycosylated systems are found to exhibit predominantly funnel-like landscapes, with a limited number of metastable states in the vicinity of the native protein fold. In contrast, glycosylation enhances landscape complexity, resulting in a wide array of relaxation timescales. Strikingly, the two glycan chains affect the landscape topography in distinct ways, despite having closely matching sequences. Dimerization couples the glycan chain dynamics, with transitions between key metastable states involving coordinated motions of both the chains. Our work illustrates that interpretation in terms of disconnectivity graphs and transition networks could reveal important insights into the organization of glycoprotein energy landscapes.

All-atom simulations with explicit solvent provide the most detailed and accurate description of dynamics and mechanisms in intrinsically disordered proteins and their condensates. However, interactions involving charged residues and ions remain a persistent source of systematic error. Here we introduce an osmometry-guided optimization strategy that directly targets residue-residue, residue-ion and ion-ion interactions. Osmotic pressure provides key experimental information on molecular interactions and can be calculated directly and rapidly from simulations, enabling efficient iterative force-field optimization. The resulting parameters improve agreement of all-atom simulations with a range of experimental data: single-molecule FRET measurements for 16 monomeric intrinsically disordered regions; NMR relaxation data for a complex between an IDP and a folded protein domain; and mean FRET efficiencies and chain reconfiguration times of IDPs in biomolecular condensates of highly charged proteins. For such condensates, simulations with an osmometry-calibrated force field provide the missing link for predicting condensate dynamics across length and time scales. The presented optimization strategy is broadly extensible to other interaction classes, including those governing protein-DNA and protein-RNA assemblies.

The fidelity of molecular dynamics (MD) simulations fundamentally depends on the quality and coverage of the ab initio data used to parameterize the underlying force field, yet the role of side-chain conformational space remains insufficiently explored. In this study, we systematically investigate how comprehensive ab initio sampling of dipeptide conformations--specifically targeting side-chain degrees of freedom--impacts force field accuracy and MD simulation predictive power. We present the Protein Unit Dataset Plus (PUD+), a 40-million-conformation quantum mechanical dataset featuring unprecedented coverage of both backbone and side-chain conformational space. Machine learning force fields trained on PUD+ and integrated into AI2BMD simulations demonstrate superior energy and force prediction accuracy, capturing high-fidelity protein folding dynamics and the conformational flexibility of long-side-chain systems. Furthermore, leveraging PUD+ to reparameterize the CMAP term of the classical ff19SB force field markedly improves the description of intrinsically disordered protein (IDP) dynamics and IDP-ligand binding. Collectively, these results demonstrate that ab initio sampling of dipeptide side-chain conformations enables high-fidelity modeling of protein dynamics across both AI-driven and classical simulation paradigms.

In linear polysomes, excluded-volume interactions among ribosomes can induce dimensional reduction of mRNA. Yet linear architectures allow steric stress to relax at open ends-- limiting how strongly crowding can remodel the mRNA's structure and dynamics. Using coarse-grained molecular-dynamics simulations, we compare circular and linear polysomes over a range of ribosome densities. Circular closure selects a predominantly quasi-planar global conformational ensemble, as indicated by a shape dimensionality dshape {approx} 2 over a range of ribosome densities. Crucially, circular topology and ribosome crowding act cooperatively to suppress structural fluctuations. While closure alone or linear crowding reduces relative global size fluctuations ({Delta}Rg/Rg) only to {approx} 0.16, their combined effect drives this fluctuation down to {approx} 0.07. Within this stabilized architecture, increasing ribosome density drives a distinct in-plane reorganization: the ring becomes more isotropic, global size fluctuations are strongly suppressed, and the scaling exponent increases toward {nu} [~=] 0.74 - 0.77, consistent with two-dimensional self-avoiding walk-like value over the accessible finite-size window, 1000 [&le;] N [&le;] 4969. Closure shortens the radius-of-gyration decorrelation time of circular polysomes by 40-fold relative to matched linear systems, reflecting the topological elimination of free ends. Within this closureselected ensemble, ribosome crowding further reduces the decorrelation time by up to 20% at the highest density. A fluctuation-informed crossover model links the density dependence of the global scaling exponent to inter-ribosomal subchain statistics. These results distinguish the geometric role of circular closure from the density-dependent steric response that it enables, revealing a confined yet dynamically responsive conformational regime for circular polysomes.